clear all; close all; clc

H = 19100e3;
Rz = 6400e3;
Tobr = 12 * 60 * 60;
Tmod = Tobr / 2;
wobr = 2*pi / Tobr;
c_light = 3e8;
f_0 = 1602e6;

dt = 5 * 60;
t = 0:dt:Tmod;

rr = [0; Rz];

D = nan(1, length(t));

rsv = (H + Rz) * ...
       [cos(wobr*t);
        sin(wobr*t)];
rv = rsv - repmat(rr, 1, length(t));
D = sqrt(rv(1, :).^2 + rv(2, :).^2);
ind = find(rv(2, :) > 0);
V = diff(D)/dt;
df = - V/c_light * f_0;

figure;
subplot(2, 1, 1)
plot(t(2:end), df, t(ind), df(ind-1));
xlabel('t, sec')
ylabel('\Delta f, Hz')

subplot(2, 1, 2)
hist(df(ind-1), 10)
xlabel('Delta f, Hz')


